// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPLINE_FITTING_H
#define EIGEN_SPLINE_FITTING_H

#include <algorithm>
#include <functional>
#include <numeric>
#include <vector>

#include "SplineFwd.h"

#include "../../../../Eigen/LU"
#include "../../../../Eigen/QR"

namespace Eigen {
/**
   * \brief Computes knot averages.
   * \ingroup Splines_Module
   *
   * The knots are computed as
   * \f{align*}
   *  u_0 & = \hdots = u_p = 0 \\
   *  u_{m-p} & = \hdots = u_{m} = 1 \\
   *  u_{j+p} & = \frac{1}{p}\sum_{i=j}^{j+p-1}\bar{u}_i \quad\quad j=1,\hdots,n-p
   * \f}
   * where \f$p\f$ is the degree and \f$m+1\f$ the number knots
   * of the desired interpolating spline.
   *
   * \param[in] parameters The input parameters. During interpolation one for each data point.
   * \param[in] degree The spline degree which is used during the interpolation.
   * \param[out] knots The output knot vector.
   *
   * \sa Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data
   **/
template <typename KnotVectorType> void KnotAveraging(const KnotVectorType& parameters, DenseIndex degree, KnotVectorType& knots)
{
    knots.resize(parameters.size() + degree + 1);

    for (DenseIndex j = 1; j < parameters.size() - degree; ++j) knots(j + degree) = parameters.segment(j, degree).mean();

    knots.segment(0, degree + 1) = KnotVectorType::Zero(degree + 1);
    knots.segment(knots.size() - degree - 1, degree + 1) = KnotVectorType::Ones(degree + 1);
}

/**
   * \brief Computes knot averages when derivative constraints are present.
   * Note that this is a technical interpretation of the referenced article
   * since the algorithm contained therein is incorrect as written.
   * \ingroup Splines_Module
   *
   * \param[in] parameters The parameters at which the interpolation B-Spline
   *            will intersect the given interpolation points. The parameters
   *            are assumed to be a non-decreasing sequence.
   * \param[in] degree The degree of the interpolating B-Spline. This must be
   *            greater than zero.
   * \param[in] derivativeIndices The indices corresponding to parameters at
   *            which there are derivative constraints. The indices are assumed
   *            to be a non-decreasing sequence.
   * \param[out] knots The calculated knot vector. These will be returned as a
   *             non-decreasing sequence
   *
   * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
   * Curve interpolation with directional constraints for engineering design. 
   * Engineering with Computers
   **/
template <typename KnotVectorType, typename ParameterVectorType, typename IndexArray>
void KnotAveragingWithDerivatives(const ParameterVectorType& parameters, const unsigned int degree, const IndexArray& derivativeIndices, KnotVectorType& knots)
{
    typedef typename ParameterVectorType::Scalar Scalar;

    DenseIndex numParameters = parameters.size();
    DenseIndex numDerivatives = derivativeIndices.size();

    if (numDerivatives < 1)
    {
        KnotAveraging(parameters, degree, knots);
        return;
    }

    DenseIndex startIndex;
    DenseIndex endIndex;

    DenseIndex numInternalDerivatives = numDerivatives;

    if (derivativeIndices[0] == 0)
    {
        startIndex = 0;
        --numInternalDerivatives;
    }
    else
    {
        startIndex = 1;
    }
    if (derivativeIndices[numDerivatives - 1] == numParameters - 1)
    {
        endIndex = numParameters - degree;
        --numInternalDerivatives;
    }
    else
    {
        endIndex = numParameters - degree - 1;
    }

    // There are (endIndex - startIndex + 1) knots obtained from the averaging
    // and 2 for the first and last parameters.
    DenseIndex numAverageKnots = endIndex - startIndex + 3;
    KnotVectorType averageKnots(numAverageKnots);
    averageKnots[0] = parameters[0];

    int newKnotIndex = 0;
    for (DenseIndex i = startIndex; i <= endIndex; ++i) averageKnots[++newKnotIndex] = parameters.segment(i, degree).mean();
    averageKnots[++newKnotIndex] = parameters[numParameters - 1];

    newKnotIndex = -1;

    ParameterVectorType temporaryParameters(numParameters + 1);
    KnotVectorType derivativeKnots(numInternalDerivatives);
    for (DenseIndex i = 0; i < numAverageKnots - 1; ++i)
    {
        temporaryParameters[0] = averageKnots[i];
        ParameterVectorType parameterIndices(numParameters);
        int temporaryParameterIndex = 1;
        for (DenseIndex j = 0; j < numParameters; ++j)
        {
            Scalar parameter = parameters[j];
            if (parameter >= averageKnots[i] && parameter < averageKnots[i + 1])
            {
                parameterIndices[temporaryParameterIndex] = j;
                temporaryParameters[temporaryParameterIndex++] = parameter;
            }
        }
        temporaryParameters[temporaryParameterIndex] = averageKnots[i + 1];

        for (int j = 0; j <= temporaryParameterIndex - 2; ++j)
        {
            for (DenseIndex k = 0; k < derivativeIndices.size(); ++k)
            {
                if (parameterIndices[j + 1] == derivativeIndices[k] && parameterIndices[j + 1] != 0 && parameterIndices[j + 1] != numParameters - 1)
                {
                    derivativeKnots[++newKnotIndex] = temporaryParameters.segment(j, 3).mean();
                    break;
                }
            }
        }
    }

    KnotVectorType temporaryKnots(averageKnots.size() + derivativeKnots.size());

    std::merge(averageKnots.data(),
               averageKnots.data() + averageKnots.size(),
               derivativeKnots.data(),
               derivativeKnots.data() + derivativeKnots.size(),
               temporaryKnots.data());

    // Number of knots (one for each point and derivative) plus spline order.
    DenseIndex numKnots = numParameters + numDerivatives + degree + 1;
    knots.resize(numKnots);

    knots.head(degree).fill(temporaryKnots[0]);
    knots.tail(degree).fill(temporaryKnots.template tail<1>()[0]);
    knots.segment(degree, temporaryKnots.size()) = temporaryKnots;
}

/**
   * \brief Computes chord length parameters which are required for spline interpolation.
   * \ingroup Splines_Module
   *
   * \param[in] pts The data points to which a spline should be fit.
   * \param[out] chord_lengths The resulting chord length vector.
   *
   * \sa Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data
   **/
template <typename PointArrayType, typename KnotVectorType> void ChordLengths(const PointArrayType& pts, KnotVectorType& chord_lengths)
{
    typedef typename KnotVectorType::Scalar Scalar;

    const DenseIndex n = pts.cols();

    // 1. compute the column-wise norms
    chord_lengths.resize(pts.cols());
    chord_lengths[0] = 0;
    chord_lengths.rightCols(n - 1) = (pts.array().leftCols(n - 1) - pts.array().rightCols(n - 1)).matrix().colwise().norm();

    // 2. compute the partial sums
    std::partial_sum(chord_lengths.data(), chord_lengths.data() + n, chord_lengths.data());

    // 3. normalize the data
    chord_lengths /= chord_lengths(n - 1);
    chord_lengths(n - 1) = Scalar(1);
}

/**
   * \brief Spline fitting methods.
   * \ingroup Splines_Module
   **/
template <typename SplineType> struct SplineFitting
{
    typedef typename SplineType::KnotVectorType KnotVectorType;
    typedef typename SplineType::ParameterVectorType ParameterVectorType;

    /**
     * \brief Fits an interpolating Spline to the given data points.
     *
     * \param pts The points for which an interpolating spline will be computed.
     * \param degree The degree of the interpolating spline.
     *
     * \returns A spline interpolating the initially provided points.
     **/
    template <typename PointArrayType> static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree);

    /**
     * \brief Fits an interpolating Spline to the given data points.
     *
     * \param pts The points for which an interpolating spline will be computed.
     * \param degree The degree of the interpolating spline.
     * \param knot_parameters The knot parameters for the interpolation.
     *
     * \returns A spline interpolating the initially provided points.
     **/
    template <typename PointArrayType> static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters);

    /**
     * \brief Fits an interpolating spline to the given data points and
     * derivatives.
     * 
     * \param points The points for which an interpolating spline will be computed.
     * \param derivatives The desired derivatives of the interpolating spline at interpolation
     *                    points.
     * \param derivativeIndices An array indicating which point each derivative belongs to. This
     *                          must be the same size as @a derivatives.
     * \param degree The degree of the interpolating spline.
     *
     * \returns A spline interpolating @a points with @a derivatives at those points.
     *
     * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
     * Curve interpolation with directional constraints for engineering design. 
     * Engineering with Computers
     **/
    template <typename PointArrayType, typename IndexArray>
    static SplineType
    InterpolateWithDerivatives(const PointArrayType& points, const PointArrayType& derivatives, const IndexArray& derivativeIndices, const unsigned int degree);

    /**
     * \brief Fits an interpolating spline to the given data points and derivatives.
     * 
     * \param points The points for which an interpolating spline will be computed.
     * \param derivatives The desired derivatives of the interpolating spline at interpolation points.
     * \param derivativeIndices An array indicating which point each derivative belongs to. This
     *                          must be the same size as @a derivatives.
     * \param degree The degree of the interpolating spline.
     * \param parameters The parameters corresponding to the interpolation points.
     *
     * \returns A spline interpolating @a points with @a derivatives at those points.
     *
     * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008.
     * Curve interpolation with directional constraints for engineering design. 
     * Engineering with Computers
     */
    template <typename PointArrayType, typename IndexArray>
    static SplineType InterpolateWithDerivatives(const PointArrayType& points,
                                                 const PointArrayType& derivatives,
                                                 const IndexArray& derivativeIndices,
                                                 const unsigned int degree,
                                                 const ParameterVectorType& parameters);
};

template <typename SplineType>
template <typename PointArrayType>
SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters)
{
    typedef typename SplineType::KnotVectorType::Scalar Scalar;
    typedef typename SplineType::ControlPointVectorType ControlPointVectorType;

    typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;

    KnotVectorType knots;
    KnotAveraging(knot_parameters, degree, knots);

    DenseIndex n = pts.cols();
    MatrixType A = MatrixType::Zero(n, n);
    for (DenseIndex i = 1; i < n - 1; ++i)
    {
        const DenseIndex span = SplineType::Span(knot_parameters[i], degree, knots);

        // The segment call should somehow be told the spline order at compile time.
        A.row(i).segment(span - degree, degree + 1) = SplineType::BasisFunctions(knot_parameters[i], degree, knots);
    }
    A(0, 0) = 1.0;
    A(n - 1, n - 1) = 1.0;

    HouseholderQR<MatrixType> qr(A);

    // Here, we are creating a temporary due to an Eigen issue.
    ControlPointVectorType ctrls = qr.solve(MatrixType(pts.transpose())).transpose();

    return SplineType(knots, ctrls);
}

template <typename SplineType>
template <typename PointArrayType>
SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree)
{
    KnotVectorType chord_lengths;  // knot parameters
    ChordLengths(pts, chord_lengths);
    return Interpolate(pts, degree, chord_lengths);
}

template <typename SplineType>
template <typename PointArrayType, typename IndexArray>
SplineType SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points,
                                                                 const PointArrayType& derivatives,
                                                                 const IndexArray& derivativeIndices,
                                                                 const unsigned int degree,
                                                                 const ParameterVectorType& parameters)
{
    typedef typename SplineType::KnotVectorType::Scalar Scalar;
    typedef typename SplineType::ControlPointVectorType ControlPointVectorType;

    typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType;

    const DenseIndex n = points.cols() + derivatives.cols();

    KnotVectorType knots;

    KnotAveragingWithDerivatives(parameters, degree, derivativeIndices, knots);

    // fill matrix
    MatrixType A = MatrixType::Zero(n, n);

    // Use these dimensions for quicker populating, then transpose for solving.
    MatrixType b(points.rows(), n);

    DenseIndex startRow;
    DenseIndex derivativeStart;

    // End derivatives.
    if (derivativeIndices[0] == 0)
    {
        A.template block<1, 2>(1, 0) << -1, 1;

        Scalar y = (knots(degree + 1) - knots(0)) / degree;
        b.col(1) = y * derivatives.col(0);

        startRow = 2;
        derivativeStart = 1;
    }
    else
    {
        startRow = 1;
        derivativeStart = 0;
    }
    if (derivativeIndices[derivatives.cols() - 1] == points.cols() - 1)
    {
        A.template block<1, 2>(n - 2, n - 2) << -1, 1;

        Scalar y = (knots(knots.size() - 1) - knots(knots.size() - (degree + 2))) / degree;
        b.col(b.cols() - 2) = y * derivatives.col(derivatives.cols() - 1);
    }

    DenseIndex row = startRow;
    DenseIndex derivativeIndex = derivativeStart;
    for (DenseIndex i = 1; i < parameters.size() - 1; ++i)
    {
        const DenseIndex span = SplineType::Span(parameters[i], degree, knots);

        if (derivativeIndex < derivativeIndices.size() && derivativeIndices[derivativeIndex] == i)
        {
            A.block(row, span - degree, 2, degree + 1) = SplineType::BasisFunctionDerivatives(parameters[i], 1, degree, knots);

            b.col(row++) = points.col(i);
            b.col(row++) = derivatives.col(derivativeIndex++);
        }
        else
        {
            A.row(row).segment(span - degree, degree + 1) = SplineType::BasisFunctions(parameters[i], degree, knots);
            b.col(row++) = points.col(i);
        }
    }
    b.col(0) = points.col(0);
    b.col(b.cols() - 1) = points.col(points.cols() - 1);
    A(0, 0) = 1;
    A(n - 1, n - 1) = 1;

    // Solve
    FullPivLU<MatrixType> lu(A);
    ControlPointVectorType controlPoints = lu.solve(MatrixType(b.transpose())).transpose();

    SplineType spline(knots, controlPoints);

    return spline;
}

template <typename SplineType>
template <typename PointArrayType, typename IndexArray>
SplineType SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points,
                                                                 const PointArrayType& derivatives,
                                                                 const IndexArray& derivativeIndices,
                                                                 const unsigned int degree)
{
    ParameterVectorType parameters;
    ChordLengths(points, parameters);
    return InterpolateWithDerivatives(points, derivatives, derivativeIndices, degree, parameters);
}
}  // namespace Eigen

#endif  // EIGEN_SPLINE_FITTING_H
